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We develop the excitation operator method, which is designed to solve the Heisenberg equation 
of motion by constructing the excitation operators. We use it to study the spin dynamics in the 
one-dimensional XXZ model. We find the difi'usive spin transport in the gapped phase at the high 
temperature limit. 
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I. INTRODUCTION 

Recently, the nonequilibrium dynamics of the zero 
and one-dimensional strongly-correlated systems attracts 
much attention due to the requirements of explaining a 
wide range of experiments in a variety of real materials. 
In spite of the intense efforts, the nonequilibrium OD and 
ID systems still reject a full understanding, due to the 
lack of a reliable numerical or analytical tool of solving 
the Schrodinger equation beyond the perturbation the- 
ory. 

The attempts of tackling this problem result in a lot of 
new nonequilibrium methods. Most of them are based on 
the developments of their corresponding versions for the 
equilibrium systems. The examples include the real time 
Quantum Monte Carlo [ll-Q], the time-dependent numer- 
ical renormalization group 0-0, M, the scattering state 
numerical renormalization group [8|, the adaptive time- 

tendent density matrix renormalization group [lol - 
the scattering Bethe-ansatz [Toj . the method based 
on integrability pol [2ll |. the real-time renormalization 
group |22h 26 | . the real time functional renormalization 
group [271 [28| and the nonequilibrium flow equation [2§- 
|33|. However, no unique method is available which can 
both give the results in a high precision and incorporate 
the physics of the problem throughout the whole param- 
eter space. It is still necessary to study new methods. 

A usual way of driving an equilibrium system out of 
equilibrium is by a quantum quench, that is to switch 
on a Hamiltonian H which is not commutative with the 
density matrix of the system at time t = 0. The dy- 
namics of a physical quantity is studied by calculating 
the expectation value of the corresponding operator with 
respect to the solution of the Schrodinger equation. Re- 
cently, the author [34[ suggested a different way of study- 
ing the nonequilibrium dynamics by solving the Heisen- 
berg equation of motion with the excitation operators. 
By decomposing the observable operator into the linear 
combination of a series of excitation operators, we cir- 
cumvent the expansion of the S-matrix in solving the 
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Heisenberg equation. This method can be applied in the 
strongly-correlated systems beyond the perturbation the- 
ory. 

In this paper, we develop a modified version of the ex- 
citation operator method and use it to study the spin dy- 
namics in the one-dimensional XXZ model, which is an 
important paradigm of the low dimensional magnetism. 
The spin dynamics in a ID system has been a field of 
active study for a long time [35l - l49| . There are still sev- 
eral open questions remaining to be settled. One of them 
is whether an integrable system in the gapped phase ex- 
hibits a ballistic spin transport at all temperatures [ssj . 
We contribute to this question by giving a negative an- 
swer. We find that the spin transport in the gapped 
phase is diffusive at the high temperature limit, contrary 
to the gapless phase. 

The plan of the paper is the following. In Sec. II we 
introduce the excitation operator method. In Sec. HI 
we define an orthonormal basis of the operator space in 
the spin-1/2 systems. The modified excitation operator 
method will be shown in Sec. IV, and the results of the 
spin dynamics in the XXZ chain will be discussed in 
Sec. V. We comment on the method in Sec. VI. 



II. SOLVE THE HEISENBERG EQUATION OF 
MOTION BY CONSTRUCTING EXCITATION 
OPERATORS 

Up to now, the analytical and numerical methods for 
studying the strongly-correlated system focus on con- 
structing the quantum state, including the ground state, 
the density matrix at finite temperature and the time- 
dependent state in nonequilibrium. They reach to the 
target by either assuming an ansatz of the needed state 
or approaching gradually to it in a process of thinning out 
of degrees of freedom in the Hamiltonian. In these meth- 
ods the fact is neglected that the quantum state contains 
all the information of the system, much more than what 
we need to calculate a physical quantity, especially a local 
one. Therefore, a question naturally arises; is it easier to 
calculate a local physical quantity than to construct the 
full state of the system? 

In this paper we study the evolution of a physical ob- 
servable in a strongly-correlated system driven out of 
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equilibrium by solving the Heisenberg equation of mo- 
tion: 



dd{t) 
dt 



i[H,d] 



(1) 



where O is the observable operator. We solve this equa- 
tion by means of constructing the excitation operators. 
The excitation operator of a Hamiltonian H is an oper- 
ator A satisfying 



[H, A] = XA, 



(2) 



where the real number A denotes the excitation energy. 
To solve the Heisenberg equation of motion, we first de- 
compose the observable operator into the linear combi- 
nation of a series of excitation operators. Such a decom- 
position has been proved to always exist (see Ref. [Ill 
for more details). Since the evolution of the excitation 
operator is trivial as A{t) = e^'^^A, we could then get the 
evolution of the observable operator. At last we calculate 
the expectation value of the operator 0{t) with respect 
to the initial quantum state at time t — Q. To do so, we 
must be aware of the initial state beforehand, which is 
often set to be an equilibrium state and can be got by 
the different numerical or analytical methods designed 
for equilibrium systems. 

We obtain the excitation operators by the method of 
undetermined coefficients. We decompose the excitation 
operators into the linear combination of a set of oper- 
ators, which form the basis of the operator space. We 

denote the basis operators as ^Oi , O2 , • ' ' ^ Then we 
calculate the commutator between the Hamiltonian and 
these basis operators. The result is generally expressed 
in a matrix form as 



H,0,\ 



(3) 



Now we suppose that an excitation operator Ai can be 
expressed as 



Ai — AijOj 



(4) 



Substituting Eqs.[3]and|3]into Eq. [21 we obtain 

^^HijAkj = XkAk,i. (5) 



These operators are then included in the basis operators. 
We repeatedly calculate the commutators between the 
Hamiltonian and the operators in the basis, until enough 
number of basis operators are obtained. 

However, in the basis obtained by this way one cannot 
guarantee that the matrix H will be a real symmetric ma- 
trix, and then does not know if it is diagonalizable. This 
problem can be solved by constructing a set of " orthogo- 
nal" basis operators. Like what we do in a vector space, 
we could define an inner product between two operators, 
written as {Oi, Oj). The inner product should satisfy all 
the properties of that defined between two vectors. Ad- 
ditionally, we set it to be always a real number. A set of 
basis operators are orthogonal to each other if and only 
if they satisfy 



(6) 



By using this relation and Eq.[3l we find that the elements 
of H can be expressed as 



(7) 



T-Lj^i, we need choose a set of ba- 



To guarantee Hi^j 
sis operators and an appropriate definition of the inner 
product, so that 



(0„[ff,6j]) = (Oj,[i7,0,]). 



(8) 



This is equivalent to say that the superoperator [H , ■] is 
self-adjoint. As will be shown next, such a set of basis 
operators can always be found in a spin-1/2 system. 

Now since H is a real symmetric matrix and its eigen- 
vectors are just the rows of the orthogonal matrix A^ 
the basis operator Oi can be decomposed into the linear 
combination of the excitation operators: 



di — Aj.iAj 



(9) 



Its expectation value at arbitrary time can be expressed 
as 



(10) 



where (Oi) is the expectation value of the basis operator 
with respect to the initial state. 



So the coefficients of A is an eigenvector of the matrix 
H with the eigenvalue Afc denoting the excitation energy. 
These coefficients can be obtained by an arbitrary diag- 
onalization algorithm. 

The critical step of our approach is the choice of the ba- 
sis of the operator space. Since the operator space should 
whatever contain the observable operator, a natural way 
to choose the basis is to take the O itself as the first 
element of the basis. We then calculate the commuta- 
tor [H , O], which will generate new operators beyond O. 



III. AN ORTHONORMAL BASIS OF THE 
OPERATOR SPACE IN SPIN-1/2 SYSTEMS 

An orthonormal basis satisfying Eq. [5] is found in the 
spin-1/2 systems, where the spins are located at each 
site on a discrete lattice. The dimension of the operator 
space at each site is four. A natural choice of the onsite 
basis operators should include the three generators of the 
SU{2) algebra and the identity operator. Additionally, 
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we must guarantee the matrix 'H to be a real one. So we 
employ the next three operators; 



1 

1 



-1 

1 



,0" 



1 

-1 



, (11) 



instead of the traditional SU{2) generators. The com- 
mutation relations between them are [6^,6^] — 26^, 
[6^,0^] = 26^ and [6^,6^] = —26^. An arbitrary basis 
operator is expressed as 



= ]^(g)6„ 



(12) 



where i denotes the lattice site and 6i takes a value in 
the set {6^,6^',o^l}. 

The inner product of two basis operators is defined as 
the trace of their product, i.e., 



Tr 



2^ 



(13) 



where L is the total number of sites in the model. It 
is not difficult to prove that such a definition satisfies 
both the positive definiteness of the inner product and 
the Eq. |8]as the Hamiltonian is self-adjoint. 

In this paper we focus on the one-dimensional XXZ 
model, because it can be strictly diagonalized in the non- 
interacting case so that we could compare the result of 
our method with the exact result. The Hamiltonian of 
the ID XXZ model is 



H 



L-l 

E 



QX QX 



AS!Sf_ 



(14) 



where the interacting strength A controls the magnetic 
order of the system at the ground state. We are inter- 
ested in the spin transport in this model. We study the 
evolution of the z-component of a spin located at site m 
denoted by {Sj^{t)). For simplicity, the initial state is set 
to be at the high temperature limit with one spin aligned 
along the positive z-axis. This spin is located at site m', 
and the \m — m'\ is the distance between the spin that we 
measure and the spin source. At the high temperature 
limit, the z-components of the spins are all zero except 
that (S^,) = ^, and all the high order spin correlation 
functions vanish. Of course, such an initial state is much 
simpler than the ground state of the system. But this 
does not mean that there is any difference when apply- 
ing our method to the zero-temperature system. Since 
in our method we build a pure operator identity, the 
additional effort that one needs at zero temperature is 
to evaluate the expectation values of the basis operators 
with respect to the ground state by employing, e.g., the 
DMRG algorithm or the Bethe Ansatz. 

By using the basis operators we re-express the Hamil- 
tonian as 



(15) 
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FIG. 1: The spin dynamics at the site m — m' as A — 0. The 
different max level Lm is adopted in the calculation. The solid 
line denotes the exact result obtained by strictly diagonalizing 
the Hamiltonian. 



The observable operator o^j is called at the zeroth level. 
We calculate the commutator between the Hamiltonian 
and it and obtain 



+ 1 r, "r?x + 1 



(16) 



The four basis operators o^jO^+i, o^oJl^+i, dl^_^dl^ and 
Om-i^m generated above are called the operators of 
the first level. We then calculate the commutators be- 
tween the Hamiltonian and each of them, which will 
generate the basis operators at the second level, and 
so on. We use the L,„ to denote the max level that 
we reach for the decomposition of the 6^. The num- 
ber of basis operators involved for the first several Lm is 
1, 5, 21, 71, 233, 729, 2223, 6573, • • • . 

Fig. 1 shows the {Sm{t)) calculated as — 4,5,6,7, 
compared with the exact result, when the interacting 
strength A and the distance \m — m'\ are both zero. The 
transient behavior of the z-component of the target spin 
is well predicted by decomposing it into the excitation 
operators. However, a strong deviation from the exact 
result is observed in the long-period behavior at a finite 
Lm. Increasing the Lm will improve the performance at 
the long period, but will also increase the difficulty in 
diagonalizing the matrix %. For example, as Lm = 7 the 
dimension of the % will be 6573! The longer the time is, 
the larger the dimension of the T-L should be to obtain the 
correct result. The error of the calculation at a finite Lm 
can be estimated by comparing the result with that after 
increasing Lm by one. As shown in Fig. 1, this gives the 
error bound. 
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IV. THE MODIFIED EXCITATION OPERATOR 
METHOD 

With the generation of the basis operators, the dimen- 
sion of the operator space and the size of the matrix % 
wih increase rapidly, and wih be soon too large to the 
ability of the most powerful tools of diagonalization. We 
conquer this problem by dividing the evolution into A'' pe- 
riods, each of which has a very short time interval. The 
observable operator in the Heisenberg picture is then ex- 
pressed as 

(17) 

where r = t/N is the time interval of a period. At 
each period the excitation operator method is employed 
to calculate the evolution of the involved basis opera- 
tors e'^'^0,e~*^'^. Since the short-time evolution can be 
obtained to a high precision by the excitation operator 
method, the d^{t) at whatever time can be obtained ex- 
plicitly as T — >■ 0. By dividing the time t into A'' intervals, 
we both improve the performance of our method at the 
long time region and avoid diagonalizing a large matrix 
by keeping the Lm at each period a small number. 

However, the exponential growth of the number of in- 
volved basis operators is still kept. Since the max level 
that we reach in the whole evolution process is propor- 
tional to N, the N cannot be very large. In practice, 
the N can be greatly enhanced by neglecting the basis 
operators with a small weight. In fact, we find that as 
decomposing the Oi(T), the weights of different basis op- 
erators may vary by several magnitudes. Most of the 
involved basis operators has a very small weight. So we 
set a small number e denoting the tolerance. At each pe- 
riod, the generated basis operators with a weight smaller 
than e are neglected. By doing so, we could control the 
growth of the operator space at a low speed. And it is 
not difficult to estimate the error caused by a finite e by 
letting g — >■ 0. 

Fig. 2 shows the result of the {S^{t)) as A = 0. The 
numerical results are found to coincide well with the ex- 
act result at the time as large as 9, when the excitation 
operator method without the time division fails. The ex- 
act result is obtained by strictly diagonalizing the Hamil- 
tonian, which can be re-expressed in terms of spinless 
fermions c-^^ through the Jordon-Wigner transformation: 

1 ^"^ 

H=-J2i4ci+i+hx.). (18) 

i=l 

The evolution of the density nm{t) = cln{t)cm{t) can be 
easily worked out, so does the S^{t). In the thermody- 
namic limit as L — > 00 and m = m', the result is 

SLit) = ijo'(i), (19) 
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FIG. 2: The simulation of the spin dynamics at A = by us- 
ing the modified excitation operator method. The results are 
compared with the exact one, denoted by the solid line. The 
Lm denotes the max level as decomposing the basis operators 
in one period. 



where Jg is the Bessel function of the first kind. 

Since for different L^s the numerical results always 
converge to the exact one as the tolerance e and the pe- 
riod T both go to zero, how to choose at each period 
is totally decided by the efficiency of the programme. We 
find that keeping a smallest number of basis operators at 
each period, i.e., choosing Lm = 1, will always save the 
computation time. The reason is that the accuracy of 
the result is mostly decided by the total number of basis 
operators involved in the calculation process. And the 
time cost in diagonalizing a matrix grows faster than the 
dimension of the matrix. So diagonalizing a lot of small 
matrices will be more efficient than diagonalizing a few 
of big ones. 

One should notice that taking a large tolerance at the 
same time as taking a too small r will cause a ridicu- 
lous result due to the Landau- Zener effect. Because the 
operator keeps almost the same when the evolution time 
T is very small, and a correspondingly large e will filter 
out all the other basis operators except itself and then 
the total operator space never grows. So as keeping the 
tolerance invariant, one should begin with a correspond- 
ingly large r. As t goes to zero, the result at the long 
time region will first converge towards the exact result, 
and then, after the r reaching the breakpoint, diverge as 
the T decreasing more. 

In the modified excitation operator method, the opera- 
tor space grows exponentially with the number of periods 
N (see Fig. 3) as the r is fixed, even after we set a non- 
zero tolerance. So the computational resource increases 
exponentially with the time, since N increases linearly 
with the time. But if the time t is fixed, increasing the 
N will not significantly increase the computation time, 
because the number of involved basis operators will de- 
crease as the time interval t decreasing. In this case, the 
tolerance e is critical to the computation time. 
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FIG. 3: The growing of the number of basis operators involved 
at each period with the period. Each data is titled by the pair 
(e,T). The interacting strength is A = 0.5. 
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FIG. 4: The spin dynamics at the site m = m' as A = 0.5, 
calculated by the modified excitation operator method and 

the excitation operator method. In the calculation using the 
modified method, the time interval at each period r is varied, 
while the tolerance e is fixed to be 10^'^. The curves titled 
Lm = 5 and Lm = 6 are the results got by the excitation 
operator method without the time division. 



V. THE DIFFUSIVE AND BALLISTIC 
TRANSPORT IN THE XXZ MODEL 

The XXZ model is integrable. As A ^ 0, the ground 
state of the Hamiltonian has been well studied by using 
the Bethe Ansatz and field theoretical methods. In the 
case of the antiferromagnetic coupling, the system expe- 
riences a quantum phase transition at the point A = 1 
from a gapless XY phase to a gapped antiferromagnetic 
phase. However, there is no method of diagonalizing the 
Hamiltonian directly. To study the non-equilibrium dy- 
namics of this system, one has to appeal to the numerical 
calculations. 

Fig. 4 shows the S^^{t) as to = to' and A = 0.5. The 
convergence of the result as r decreasing has been ob- 



FIG. 5: The comparison of the spin dynamics at the site 
m = m' at different A. In the calculation, the tolerance is set 
to be e = 10~^ and the time interval is set to be t = 0.04. 
An intermediate quasi-steady regime is found at A = 1.5. 
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FIG. 6: The variance functions at different A, plotted in the 

logscalo. The variance function is linear in the gapped phase 
as A = 1.5, while quadratic in the gapless phase. 



served as the tolerance e = 10""^. We also compare the 
result with that without the time division, in which a 
good estimation of the error can be obtained. We con- 
clude that the result as e — 10^ ^ has been very close to 
the real one. Decreasing e more will not result in a signif- 
icant enhancement to the accuracy. And this conclusion 
is not changed as the distance |to — to'| increasing. 

We study the spin relaxation at different A (see Fig. 5). 
An obvious difi'erence is observed between the gapless 
XY phase and the gapped anti-ferromagnetic phase: in 
the gapped phase there is an intermediate quasi-steady 
regime, which is not seen in the gapless phase. The quasi- 
steady S^, at the intermediate time before decaying to 
zero at the long time limit can be understood as the 
result of the antiferromagnetic ordering. A short-time 
local antiferromagnetic ordering around to' is induced by 
the initial spin, which again stabilize it. 

The persistence of the spin in the gapped phase indi- 
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cates a shift from the ballistic spin transport in the gap- 
less phase to the diffusive transport. For quantitatively 
distinguishing the ballistic and diffusive spin transport, 
we turn to the time-dependent variance of the spin den- 
sity, defined as 

<^''{t)^Y.{m^^if{nm{t)-n„,{0)), (20) 

following the Ref. Here n,„(t) = {S^{t)) + i is the 
particle density function, and /i the first moment of it. In 
the thermodynamic limit, there is always fj, — m' . Espe- 
cially, the variance function at A = is strictly calculated 
to be o'^(i) = X]m=i "^^"^mlO- general, the variance 
will be a linear function of the time in a diffusive trans- 
port, while a quadratic function in a ballistic transport. 
We compare the variance functions at different A. We 
observe the obvious character of the diffusive transport 
in the variance function at A = 1.5. At the same time, 
the variance function at A = 0.5 looks similar to the 
noninteracting one. 



its focusing on solving the Heisenberg equation of mo- 
tion instead of the Schrodinger equation. In this aspect 
it is similar to the real time flow equation approach. It 
will result in a pure operator identity. An extra equi- 
librium method of calculating the initial state should be 
employed if the initial state is not trivial. The modified 
excitation operator method is suitable for a local observ- 
able, e.g., a local spin, because it save the time of getting 
the knowledge of the full quantum state in a high dimen- 
sional Hilbert space. It gives a very accurate result in an 
extremely short time in the transient regime. However, 
the computational resource grows exponentially as the 
time increasing. So it is difficult to calculate the physical 
quantity at the long time period by using this method. 
Addtionally, it is natural to perform this method in the 
thermodynamic limit, because the basis operators gener- 
ated in course of time are all the local operators. Then 
the size of the system is not prerequisite to this method. 



VI. COMMENTS 



The modified excitation operator method is distin- 
guished from most of the other real time methods by 
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